%
%Calculate heartrate
%
function [HR,Real_HR] = hrc(Record,thres)
if nargin == 1
    thres = 1650;
end
    %% Setting threshold
    %Max_HR = Maximum recorded data magnitude
    %Min_HR = Minimum recorded data magnitude
    Max_HR = max(Record(thres:length(Record)-thres));
    Min_HR = min(Record(thres:length(Record)-thres));
    
    %% Threshold calculation
    %Thres_HR_up = upper threshold setting
    %Thres_HR_low = lower threshold setting
    Thres_HR_up = Max_HR - ((Max_HR-Min_HR) * 0.2);
    Thres_HR_low = Max_HR - ((Max_HR-Min_HR) * 0.8);
    
    %% Heartrate algorithms
    %
    % By setting the upper and lower threshold, the counting of heartbeat 
    % can be calculated by a 2-stage process
    %
    % Stage 1: Before > Upper threshold
    % Stage 2: After < Lower threshold
    %
    % Two stages will repeat again until all the samples have been
    % tranformed and calculated. Two-threshold approach is implemented for
    % erroneous calculation of heartrate with small fluctuation of signals
    %
    % Parameters:
    % switch_1 and index_1: temporary parameter used for
    %                       identification and indexing respectively
    %
    % time_stored: a matrix to stored the converted each magnitude of
    %              signal in every sampling period
    % 
    % time_temp_int: a temporary variable to store the "initial" position
    %                of each cycle
    %
    % time_temp_fin: a temporary variable to store the "final" position of
    %                each cycle
    %
    % switch_1 = 1: for initial recording
    
    switch_1 = 1;
    index_1 = 1;
    time_stored = zeros(1,1);
    time_temp_int = 0;

for read_HR = 1:1:length(Record)
    %% Stage 1
    %
    % Stage 1 has been defined in two conditions
    %
    % Condition 1: The first cycle of calculation
    %   Because the computer would not know the first cycle starting point
    %   is a curve which is increasing or decreasing, we have to defined
    %   the first counting would start from the upper threshold.
    %   
    %   The computer would only record position when switch_1 = 1
    %
    %   On the other hand, time_temp_int should be an empty variable before
    %   entering the loop. We have to set time_temp_int in some values
    %   before the calculation loop starts
    %
    %   switch_1 = 0 refers to the switch telling the computer we have
    %   recorded the first position, any increasing value of threshold 
    %   would not record again
    %
    
    if (Record(read_HR) > Thres_HR_up && switch_1 == 1)
        if(read_HR == 1)
            time_temp_int = read_HR;
            switch_1 = 0;
        
        %    
        % Condition 2: After the first cycle
        %   time_temp_fin: recording the final position
        %
        %   time_temp_int is not an empty variable anymore, swapping of
        %   intial and final position should proceed for future calculation
        %   
        else
            time_temp_fin = read_HR;
            time_stored(index_1) = 1/((time_temp_fin-time_temp_int)*(1/553.6));
            switch_1 = 0;
            index_1 = index_1 + 1;
            time_temp_int = time_temp_fin;
        end
    end
    
    %% Stage 2
    %
    % If the computer has passed lower than the lower threshold, 
    % the computer would start to record another position which is just
    % same magnitude of the time_temp_int.
    %
    % switch_1 = 1 telling the computer start to find the 
    % time_temp_fin position which are required for calculation
    %
    if(Record(read_HR) < Thres_HR_low && switch_1 == 0)
        switch_1 = 1;
    end
end

  %% Real time-data processing
  % This process is to create a matrix of data recording the [time,
  % Real_HR] matrix
  %
  % time = estimated real time
  % Real_HR = By difference of real time heartrate time_stored
  %
  % HR = mode of heartrate
  % Mode is used to avoid too large transient response noise adheres which
  % affecting the calculation
  %
  Real_HR(1:length(time_stored)) = time_stored(1:length(time_stored)) * 60;
  
  %% Checking for accurate heartbeat
  % Purpose: To eliminate any vibration during measurement
  % 
  % Approach: Record a set of new data and eliminate data when the second
  %           data is three times larger than the first data
  %
  %% "Mean" heartrate calculation and selection
  % Mean calculation when standard deviation of heartrate is small
  % Mode calculation when standard deviation of heartrate is large
  
  if std(Real_HR(3:length(Real_HR))) <= 2
      HR = mean(time_stored) * 60;
  else
      HR = median(time_stored) * 60;
  end
  
  return
  

    